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SUMMARY 

We extend the Balancing Domain Decomposition by Constraints (BDDC) method to flows in porous media 
discretised by mixed-hybrid finite elements with combined mesh dimensions. Such discretisations appear 
when major geological fractures are modelled by ID or 2D elements inside three-dimensional domains. In 
this set-up, the global problem as well as the substructure problems have a symmetric saddle-point structure, 
containing a ‘penalty’ block due to the combination of meshes. We show that the problem can be reduced 
by means of iterative substructuring to an interface problem, which is symmetric and positive definite. The 
interface problem can thus be solved by conjugate gradients with the BDDC method as a preconditioner. 
A parallel implementation of this algorithm is incorporated into an existing software package for subsurface 
flow simulations. We study the performance of the iterative solver on several academic and real-world 
problems. Numerical experiments illustrate its efficiency and scalability. Preprint version. 
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1. INTRODUCTION 

A detailed description of flow in porous media is essential for building mathematical models 
with applications in, for example, water management, oil and gas recovery, carbon dioxide (CO 2 ) 
sequestration or nuclear waste disposal. In order to set up a reliable numerical model, one needs to 
have a good knowledge of the problem geometry and input parameters. For example, the flow of 
water in granite rock, which is a suitable site for nuclear waste disposal, is driven by the complex 
system of vugs, cavities and fractures with various topology and sizes. These alter the effective 
permeability, and therefore should be accurately accounted for in the numerical model. There are 
two main approaches: either the fractures are considered as free-flow regions, or the fractures 
contain debris and are also modelled as porous media with specific permeabilities. In the first case, 
a unified approach to modelling free-flow and porous media regions can be provided by the so called 
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Stokes-Brinkman equation, which reduces to either the Stokes or Darcy model in certain parameter 
limits, e.g., within the Multiscale Mixed Finite-Element (MsMFE) framework [1]. In this paper, we 
consider the latter case, and apply Darcy’s law to the flow in the reservoir and in the fractures as well; 
see [2] for a related approach. In either case, the preferential flow in large geological dislocations 
and their intersections should be considered as two- and one-dimensional flows, respectively. Due 
to the quite complex structure of the domains, the discretisation is performed using finite element 
methods (EEM). The resulting meshes are therefore unstructured, and they combine different spatial 
dimensions (line elements in ID, triangles in 2D, and tetrahedrons in 3D). The systems of linear 
equations obtained from the EEM discretisation are often very large, so that using direct methods 
is prohibitive and iterative solvers are warranted. The systems are typically also bad-conditioned 
due to the mixing of spatial dimensions, large jumps in permeability coefficients and presence of 
elements of considerably different sizes, and so they are challenging for iterative solvers as well. 

The matrices have a saddle-point structure 


A 

B 



( 1 ) 


where A is symmetric positive definite on the kernel of B, C is symmetric positive semi-definite, 

_ T _ 

and it is positive definite on the kernel of B . The ‘penalty’ block C 7 ^ 0 arises from connecting 
meshes of different spatial dimensions. The iterative solution of systems with this structure is a 
frequently studied topic; see, for example, [3, 4, 5, 6 ], the monographs [7, 8 ] or [9, Chapter 9] 
and the references therein. However, efficient methodologies for solving saddle-point problems are 
typically problem dependent. 

In this paper, we develop a robust and scalable solver for linear systems with the saddle-point 
structure as in (1) with the block C either zero or nonzero. The solver is tailored to the mixed- 
hybrid formulation of flow in porous media using the lowest order Raviart-Thomas (RTq) finite 
elements with combined mesh dimensions (ID, 2D and 3D). In particular, we adapt the Balancing 
Domain Decomposition by Constraints (BDDC) method to this type of problems. 

The BDDC method is currently one of the most popular methods of iterative substructuring. It has 
been proposed independently in [10, 11, 12]; see [13, 14] for the proof of equivalence. Even though 
BDDC has been originally formulated for elliptic problems, it has been successfully extended, 
for example, beyond elliptic cases [15, 16] and to multiple levels [17, 18]. An optimal set-up has 
been studied in [19, 20, 21, 22]. A closely related BDDC preconditioner for vector field problems 
discretised with Raviart-Thomas finite elements has been studied in [23]. 

We are interested in applications of the BDDC method to saddle-point problems. If C = 0 in (1), 
one possible approach is to use an algebraic trick and constrain the iterative solution of the indefinite 
problem into a balanced subspace, which is sometimes also called benign, where the operator 
is positive definite; see [15] for the Stokes problem, and [5, 24, 25] for flow in porous media. 
However, due to the mixed-hybrid formulation and possible coupling of meshes with different 
spatial dimensions, C 7 ^ 0 in general, and we wifi favour an alternative, dual approach here. 

Our methodology is as follows. The mixed-hybrid formulation [26, 27] is used in order to modify 
the saddle-point problem to one which is symmetric and positive definite by means of iterative 
substructuring. In particular, we introduce a symmetric positive definite Schur complement with 
respect to interface Eagrange multipliers, corresponding to a part of block C. The reduced system is 
solved by the preconditioned conjugate gradient (PCG) method, and the BDDC method is used as 
a preconditioner. From this perspective, our work can be viewed as a further extension of [ 6 ]. Our 
main effort here is in accommodating the BDDC solver to flows in porous media with combined 
mesh dimensions. In addition, the presentation of the BDDC algorithm is driven more by an efficient 
implementation, while it is more oriented towards underlying theory in [ 6 ]. We take advantage 
of the special structure of the blocks in matrix (1) studied in detail in [26, 28, 29]. In particular, 
the nonzero structure of block C resulting from a combination of meshes with different spatial 
dimensions is considered in [30]. We describe our parallel implementation of the method and 
study its performance on several benchmark and real-world problems. Another original contribution 
of this paper is proposing a new scaling operator in the BDDC method suitable for the studied 
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problems. We note that if there is no coupling of meshes with different spatial dimensions present 
in the discretisation, the block C = 0 in (1) and our method is almost identical to the one introduced 
in [6]. 

The paper is organised as follows. In Section 2, we introduce the model problem. In Section 3, we 
describe the modelling of fractured porous media and combining meshes of different dimensions. 
In Section 4, we introduce the substructuring components and derive the interface problem. In 
Section 5, we formulate the BDDC preconditioner. In addition, the selection of interface weights for 
BDDC is studied in detail in Section 6. In Section 7, we describe our parallel implementation, and in 
Section 8 we report the numerical results and parallel performance for benchmark and engineering 
problems. Finally, Section 9 provides a summary of our work. 

Our notation does not, for simplicity, distinguish between hnite element functions and 
corresponding algebraic vectors of degrees of freedom, and between linear operators and matrices 
within a specihc basis—the meaning should be clear from the context. The transpose of a matrix 
is denoted by superscript ^ and the energy norm of a vector x is denoted by \\x \\m = y'x^Mx, 
where M is a symmetric positive definite matrix. 


2. MODEL PROBLEM 


Let O be an open bounded polyhedral domain in M^. We are interested in the solution of the 
following problem, combining Darcy’s law and the equation of continuity written as 


k ^u + Vp = —Vz in fl, (2) 

V • u = / in D, (3) 

P = PN ond^N, (4) 

un = 0 onSDs, (5) 


subject to boundary conditions on 90 = where stands for the part of the 

boundary with natural (Dirichlet) boundary condition, and d^E for the part with essential 
(Neumann) boundary condition. In applications, the variable u describes the velocity of the fluid 
and p the pressure (head) in an aquifer O, k is a symmetric positive definite tensor of the hydraulic 
conductivity, — Vz = (0,0, —1)^ is the gravity term, and n is the outer unit normal vector of 90. The 
term Vz is present due to the fact, that u satisfies u = —'kVph, where pt = p + z is fhe piezometric 
head. Lor a thorough discussion of application background we refer, e.g., to monographs [31, 32]. 

Let T be the triangulation of domain O consisting of Ne simplicial elements with characteristic 
size h. We introduce a space 

H(0; div) = {v : v G L^(0); V • v G L^{fl) and v • n = 0 on 90^;} , (6) 

equipped with the standard norm. Let V c H(0,div) be the space consisting of the lowest order 
Raviart-Thomas (i?To) functions and let Q c (O) be the space consisting of piecewise constant 
functions on the elements of the triangulation T. We refer, e.g., to monograph [33] for a detailed 
description of the mixed finite elements and the corresponding spaces. 

In the mixed finite element approximation of problem (2)-(5) we look for a pair {u,p} G V x Q 
that satisfies 


u-vdx— / p\7-vdx = — / pN'v-nds— / Vzdx, Vv G V, 

J Cl J OQiv Cl 


(7) 


— / q\7 ■ udx = — / fqdx, '^qGQ. 


( 8 ) 


In the discrete formulation, we need pjy and / only sufficiently regular so that the integrals in the 
weak formulation (7)-(8) make sense, namely p^ G f & (f^)- 

Next, we recall the mixed-hybrid formulation. It was originally motivated by an effort to modify 
the saddle-point problem (7)-(8) to one which leads to symmetric positive dehnite matrices. 
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Nevertheless, this formulation is also suitable for a combination of meshes with different spatial 
dimensions, which will be described in detail in the next section. 

Let IF denote the set of inter-element /acei' of the triangulation T- We now introduce several 
additional spaces. First, let us dehne the space by relaxing the condition of continuity of the 
normal components in the space V on inter-element boundaries F. More precisely, we dehne local 
spaces V* for each element T* s T, * = 1,..., iV^;, by 

V* = {vGH(r;div):vei?To(r)}, (9) 

and put x • • • x Next, we dehne the space of Lagrange multipliers A consisting 

of functions that take constant values on individual inter-element faces in F, 


K={\eL^{F):\ = w n|^, v G V} . 


( 10 ) 


In particular, A = 0 on 90 for any A G A. 

In the mixed-hybrid finite element approximation of problem (2)-(5), we look for a triple 
{u,p, A} G X Q X K that satishes 


Ne r 


E 


• V dx — 


UT' 


/ pV dx + / A(v • n)|ari ds 
Jt' JdTFsn 


( 11 ) 


\dn 

Ne 


Pnn 




nds — ^^ / Vzdx, Vv G V 


Ne 

-E 

2 = 1 


Ne r p 


gV ■ udx 
p(u • n)|aTi ds 


aT'\an 


= - / fqdx, Vg G Q, (12) 
Jn 

= 0, V/r G A. (13) 


Equation (13) imposes a continuity condition on the normal component of the velocity (also 
called normal flux) u • n across F which guarantees that u G V. This condition also implies the 
equivalence of the two formulations (7)-(8) and (11)-(13). We note that the Lagrange multipliers A 
can be interpreted as the approximation of the trace of p on F, see [34] for details. 

Let us now write the matrix formulation corresponding to (11)-(13) as 

'A B'^ Bl 
BOO 
_ B^ 0 0 

It is important to note that A is block diagonal with Ne blocks, corresponding to elements T*, 
i = 1,., Ne, and each of the blocks is symmetric positive dehnite, cf. the hrst term in (11). It was 
shown in [28] that the system of equations (14) can be reduced (twice) to the Schur complement 
corresponding to the Lagrange multipliers A and solved efficiently by a direct or iterative solver. 
Here, we will look for an efficient solution of a slightly modified, and in general also block dense, 
system which is introduced in the next section. 



u 


9 


p 

= 

f 


A 


0 


(14) 


3. MODELLING OE ERACTURES 

In this section, we recall the main ideas of the discrete model of the flow in fractured porous 
media that is based on connection of meshes of different dimensions as described in [30]. Eet us 
denote the full domain by H 3 = fl. Next, consider lower-dimensional domains c d = 2,3, 
such that H 2 consists of polygons and Hi consists of line segments. We will also assume that 
c dfl 2 C OHa. The first condition requires that a domain of a lower dimension cannot poke 
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out of the domain of higher dimension, while the second condition prevents domains of lower 
dimension from having boundaries in the interior of domains of higher dimension. We impose these 
conditions to avoid technical difficulties in the analysis. However, numerical evidence suggests that 
these conditions are not necessary, and in fact, they are not satisfied for the real-world problems 
presented in Section 8.2. 

For every dimension d = 1,2,3, we introduce a triangulation Td of the domain Vtd that consists 
of finite elements T^, f = 1 ,..., and satisfies the compatibility conditions 

C Td. where = |J \ d^d, (15) 

k 

T^-i n dT^ is either T^_i or 0 , (16) 

for every i £ j £ N^}, and d = 2,3. This means that elements of a lower 

dimension match faces of elements of the higher dimension. 

We consider equations (3)-(5) on the domains d = 1,2,3, completed by a slight modihcation 
of the Darcy’s law (2): 

= ( 17 ) 

Od 

where stands for the velocity integrated over the cross-section for d = 1 , 2 , i.e. the units of U 3 , 
U 2 , and ui are ms“^, m^s“^, and m^s“^, respectively. In addition, 83 = 1 , 82 is the thickness of a 
fracture, and di is the cross-section of a ID preferential channel. The effective fluid source /2 on D 2 
is given as 

/2 = <^2/2+ U;^ • n++ U3 • n", ( 18 ) 

where /2 is the density of external fluid sources, and the normal fluxes from the two faces of the 
3D continuum surrounding the fracture are given through the Robin (also called Newton) boundary 
conditions 


u+• n+= (T+(p+-P 2 ), (19) 

U 3 • n" = (T;^(P 3 -P 2 )- ( 20 ) 

In the last formula, (T 3 ^~ > 0 are the transition coefficients (cf. [ 2 ] for possible choices) and p^, p^ 
are the traces of pressure p 3 on the two sides of the fracture. The effective fluid source fi on Hi is 
similar, 

/i =^l/l+^u^n^ ( 21 ) 

k 

where fi is fhe density of external fluid sources. In the 3D ambient space, the ID channel can be 
connected to k faces of 2D fractures. Thus 

\il ■ = a^ipl - Pi) (22) 

is the normal flux from the connected fracture k, > 0 is the transition coefficient, and p^ is the 
trace of pressure p 2 on the face of fracture k. 

In the following, we describe the discrete mixed-hybrid formulation of the problem. The 
formulation and discussion of the continuous problem can be found in [30]. Let us consider spaces 

]\rd- 

i\^ 

V-I=vr'X V 2-1 X V3-\ V-1= Q = QixQ2xQ3,Qd = L^nd). (23) 

i=l 

For the definition of the space A, we cannot follow (10) directly, since e.g. on H 2 , we need to 
distinguish values of A 3 on two sides of a fracture. Thus, we introduce a separate value for every 
non-boundary side of every element: 

A(T*) = {a G L^dTl \ d^d) : A = V v G V^} , (24) 
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where is defined in the same way as the space V but on the domain Further, we identify 
values on faces/points that are not aligned to the fractures/channels: 




Aw = 


{ A e J] A(r*); XU. = on face F = dT^ dT^ if F n = 0}. (25) 


i=l 


Finally, we redefine A = Ai x A 2 x A 3 . In the mixed-hybrid finite element approximation of the 
flow in fractured porous media we seek a triple {u,p, A} S x Q x A that satisfies 


a (u, v) + 6 (p, v) + bjr (A, v) = {g, v), 
b (u, q)-c {p, q) - cjr {q, A) = (/, q) , 
(u, p) - cjr {p, p)-c (A, p) = 0, 


Vv G V-\ 
Vg G (5, 

V/r G A, 


with 


3 

a(u,v) 

d—1 i—1 

3 Ni 

Mu,g) = -EE 

i—1 

3 U 

b^ (u, A) = E E 

d—1 i—1 

3 U 

c(p,g) = EE 

d=2 i=l 
3 U 

C-T (P> M) = - E XI 

i^l 

q 

c(A,.) = EE 

d=2 i=l 
3 n; 


IT\ 


—• \ddx 
5r 


/ gw (V • Urf) dx 

dn 

Alarj (ud ■ n) ds 


idT^\dnd 


<^d Pd-i qd-i ds 


f dT^D^d-i 


' dT^nQd- 


ad Pd-iPd ds 


' dT^nQd-1 


(Td XdPdds 


3 


(g,v) = -EE/. (v-n) ds-EE / 

3 


Vz dx, 


(/, g) = - E / 


The system (26)-(28) now leads to the matrix form 


A 



B 

-C 


Bt 

-Cr 

-C 



u 


9 


P 

= 

f 


X 


0 


(26) 

(27) 

(28) 


(29) 

(30) 

(31) 

(32) 

(33) 

(34) 

(35) 

(36) 


(37) 


We note that (37) is related to (26)-(28) in the same way as (14) is related to (11)-(13). The 
main difference in the structure of the matrices between (37) and (14) is the additional block 


C = 


C U 


related to the normal fluxes between dimensions and arising from (19)-(20) 


U C 

and (22) via (32)-(34). In particular, the modified right-hand side of the continuity equation for 
two-dimensional and one-dimensional elements, /2 and /i, incorporates pressure unknowns on 2D 
and ID elements, and traces of pressure on 3D and 2D elements at the fracture, which are nothing 
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but the Lagrange multipliers on 3D and 2D elements in the mixed-hybrid method. Consequently, 
p+/- = A+/” in (19)-(20), and in (22). 

Assuming 6 ^ bounded and greater than zero, and using the fact that Ik^ corresponds to a symmetric 
positive definite matrix, we see from (29) that block A in (37) is symmetric positive dehnite. Block 
C is symmetric positive semi-dehnite since 


3 


c{p,p) + 2cjf(p, A) + c(A, A) = EE 




’ dTJ:nQd-i 


crd{Pd-i - Ad)^ ds 


(38) 


The following theorem is a standard result, e.g. [33, Theorem 1.2]. Here, we rewrite it in 
a form suitable for our setting and we verify the assumptions for the specific blocks of the 
matrix in (37). We will further denote Q = Q x A, p = {p, X) G Q, q = {q, p) £ Q, and 5(u, q) = 

b{u,q) + bjr{u,fi). 

Theorem 1 

Let natural boundary conditions (4) be prescribed at a certain part of the boundary, i.e. dTl ^4 ^ 0 
for at least one d G {1,2,3}. Then the discrete mixed-hybrid problem (37) has a unique solution. 


Proof 

Let us hrst investigate the structure of the matrix in (37) more closely. Let us number the unknowns 
within each block of (37) with respect to spatial dimension d G {1,2,3}. The matrix then takes the 
form of 9x9 blocks. 



All 



Bfi 


tdT 

■ 



A22 



B ^2 

-dT 

^T.22 





A33 


^ 3^3 


rT 

^ 4,33 


Bii 

B22 

B33 

-Cn 

— C22 

pT 

^TA2 

piT 

'■^ 4,23 


B441 









B442 


— Cjrp2 


— C22 





Sjr ,33 


— Cjr 43 


— C33 _ 

Suppose for a 

moment that we 

solve a 

problem only on 

domain d € 

1,2,3 (i.e. 


(39) 


= 0 for 

i d). If no natural boundary conditions are imposed, there is a single —1 entry on each row 
of and a single -1-1 entry on each row of Since Tld is a simply connected set, the 

_ X 

matrix = [ sL ] has a nontrivial nullspace of constant vectors. Enforcing natural boundary 
condition on a part of Qd changes the - 1-1 value on the corresponding row of matrix to 0 , in 

—T 

which case has only a trivial nullspace, i.e. full column rank (see e.g. [33, Section IV. 1] or [26, 
Lemma 3.2]). 

The nullspace becomes more complicated for domains with fractures, in which case Tld typically 
has more simply connected components separated by fractures (cf. Fig. 1). Let us denote them D):, 
k = 1,... ,nc, regardless of the dimension. In particular, k = 1,... ,nci will be components 
without natural condition boundary, i.e. dfl'j. n = 0, while for k = Ud + I, ■ ■ ■ ,nc we get 
components with prescribed natural boundary condition. We also denote Xk ^ Q the characteristic 
vector of the component that takes value - 1-1 for degrees of freedom associated with elements 
and faces of the component With such notation, the basis of the nullspace of the whole matrix 


B^ = 




Bi 


UT 


Bh 


tdT 


tdT 


(40) 


consists of characteristic vectors Xk^ ^ .. .nd and has dimension ric 
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Next, we show that matrix C is not only symmetric positive semi-dehnite, as seen from (38), but 

—T _ —T 

also positive dehnite on null B . To this end, take p G null B , a vector that is piecewise constant on 
components, having value p^. ^ on the component of dimension d for A: = 1,..., rid, and value 
Pi_ ^ = 0 for other components. Then jFCp = 0 implies p = ^ = 0. Indeed, every component 

of dimension d = 2,3 has some component fij of dimension d — 1 on its boundary, and therefore 
all p/^ ^ have the same value, cf. (38). This value is zero, since there is at least one component with 
natural boundary condition. 

Applying the congruence transformation. 


r —Ti 

A B 


I 


'A 

B -C, 


BA-^ I 


-( 


^B^ + C) 




I A-^B 
I 


(41) 


_ _ 'J' 

Matrix A is symmetric positive definite (SPD) from (29) and therefore BA~^B is SPD on range 

_ _ 'J' _ _ 'J' _ 

B, which is the orthogonal complement of the nullspace of B . Thus BA~^B + C is SPD on 
whole Q. From the Sylvester law of inertia, the number of positive, negative and zero eigenvalues is 
preserved by the congruence transformation. Since the block diagonal matrix on the right-hand side 
of (41) has only (strictly) positive and (strictly) negative eigenvalues, the matrix on the left-hand 
side also must be nonsingular, and problem (37) has a unique solution. 

□ 




Figure 1. Example of a two-dimensional problem: even single fracture gives rise to two components of the 
two-dimensional mesh (left); example of three (shaded) triangles sharing a single Lagrange multiplier in 3D 

(right). 


4. ITERATIVE SUB STRUCTURING 

Eor our purposes of combining meshes with different spatial dimensions, we dehne substructures 
as subsets of hnite elements in the mesh rather than parts of a physical domain, cf. [9]. 

To begin, let us dehne the combined triangulation 7i23 as the union of triangulations for each 
spatial dimension, i.e., 7i23 = Ti U72 UTs. The triangulation 7i23 is subsequently divided into 
substructures Ns- Note that in general a substructure can contain hnite elements of 

different dimensions. We dehne the interface T as the set of degrees of freedom shared by more than 
one substructure. Note that T c A for the current setting. Thus, let us split the Eagrange multipliers A 
into two subsets. Eirst, we denote by Ar the set of multipliers corresponding to the interface T. 
The remaining multipliers, corresponding to substructure interiors, will be denoted by A/. The 
interface T* of substructure U* is dehned as a subset of T corresponding to Next, let Ap be 
dehned as the space of Eagrange multipliers corresponding to T*, i = 1,..., Ns, and dehne a space 

Ar = A^ X ••• X A^y^. (42) 
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The substructure problems are obtained by assembling contributions of finite elements in each O*, 


A* 

RiT 

TjiT 

Bjrj 

TjiT n 


■ u* ■ 


r 5* 1 


-c* 

r<iT 

r<iT 

~ gTT 


p* 


r 


-c>./ 

r’i 

W/ 

r^iT 

— Opp 




0 

B}f,t 


/^i 

Opp 

Opp J 


L A^ J 


. 0 . 


i = l,...,Ns, (43) 


where the blocks A'^ are block diagonal with blocks corresponding to finite element matrices, and the 
blocks C‘ ^ 0 only if the substructure O* contains coupling of elements of different dimensions. We 
note that the global problem (37) could be obtained from (43) by further assembly at the interface. 

In the iterative substructuring (see e.g. [9]), we first reduce the problem to substructure interfaces. 
In our context we can eliminate normal fluxes, pressure unknowns, and Lagrange multipliers at 
interiors of substructures, and we can define the substructure Schur complements S'* : Ap i-a Ap, 
i = 1 ,..., Ns, formally as 






■ A* 

RiT 

tdIT 

Bj-J 

-1 

Bj^X 

5* = C^p + 

TDi /^i 

r^i 

Or/ 



-c* 

r-<iT 


CiT 





- B^j 


1 

1 


1 

1 

1_ 


However, in an implementation of a Krylov subspace iterative method, we only need to compute the 
matrix-vector product S*Ap for a given vector Ap. Therefore, the matrix is not constructed explicitly, 
and the multiplication is obtained as follows. 

Algorithm 2 

Given Ap € Ap, determine S*Ap S Ap in the following two steps: 


1. solve a local Dirichlet problem 


■ A* 

RiT 

jDiT "1 

Bx,i 




tdIT 


-c* 

r<iT 



= — 

r^iT 

- B^xj 

fyi 

1 

1_ 


1 

J*. 

1_ 


1 

1 

1_ 


(45) 


2. perform two sparse matrix-vector multiplications 


! ~ 



w* 

\ 

-C^pA*. + 

TDi (Di 

^TX ^xx ^r/ 


g* 


V 



/ 


(46) 


Next, let Ap be the space of global degrees of freedom, such that the values of degrees of freedom 
from Ar coincide on T. The vectors of the local substructure degrees of freedom Ap G Ap and the 
vector of the global degrees of freedom Ar G Ar are related by 

A*p = i?*Ar, i = l,...,Ns, (47) 

where ii* are the restriction operators. More specifically, each i?* is a 0-1 matrix such that every row 
contains one entry equal to one, and = I. The global Schur complement 5 : Ap —Ap can 

be obtained as 

Ns 

S = (48) 

i=l 

Equation (48) represents the formal assembly of substructure Schur complements into the global 
Schur complement. The global Schur complement system, which we would like to solve iteratively, 
reads 

S\r = b, (49) 
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where vector b = is obtained from substructure reduced right-hand sides 





■ A* 

BiT 

TDiT 1 

-1 


TDi 

r^i 

Op/ 



-C* 

r-<iT 


r 




. ^T.1 


'-'77 J 


0 


In our implementation, we factor and store the matrices from (45). The factors are then used to 
compute the vectors 6* in (50), and especially in Algorithm 2, which is used in connection to 
formula (48) to evaluate Ar —>^ S\t within each iteration of a Krylov subspace iterative method. 
This algorithm allows a straightforward parallel implementation. After an approximate solution 
of (49) is found, solution in interiors of substructures, including all primal variables, is recovered 
from (43) using the factors from (45). 

Remark 3 

There are other ways to derive the interface problem (49). The authors of [6, 34] use a mixed- 
hybrid formulation with Lagrange multipliers introduced only at the interface T as their starting 
point. While problem (49) is equivalent to the interface problems considered in [6, 34], the 
substructure problems therein have a different structure from (43). In particular, there are no 
blocks corresponding to A}, and the matrices corresponding to A* are no longer block diagonal 
element-wise. Next, the authors of [28, 30] construct explicit Schur complement with respect to the 
whole block of Lagrange multipliers A and they show that, due to the special structure of A, the 
complement is both sparse and reasonably cheap to construct. If this was performed substructure 
by substrucmre, this could be again seen as an intermediate step in obtaining problem (49) by 
additional elimination of the interior Lagrange multipliers A/. However, this would again lead to 
different substructure problems based on explicit local Schur complements. 

The following result allows an application of the BDDC method to problems with fractures. 

Theorem 4 

Let natural boundary conditions (4) be prescribed at a certain part of the boundary, i.e. dTl 7 ^ 0 
for at least one d € {1,2,3}. Then the matrix S in (49) is symmetric positive definite. 


Proof 

_ _ 'J' _ 

Using the notation of (40)-(41), let us introduce a matrix S = BA~^B -|- C. The matrix S 
is symmetric positive definite by Theorem 1 (and its proof). Applying another congruence 
transformation to S and denoting the rows corresponding to the interface Lagrange multipliers by 
subscript r and the interior by /, we obtain 


Sii 



/ 


Sii 


f s:-! CT ' 

Sti 

5rr 


SriSfl I 


Nrr — StiSjXti. 


/ 


Since the matrix on the left-hand side is SPD, both diagonal blocks Sn and Nrr — StiSJ^Syi are 
also symmetric positive definite from the Sylvester law of inertia. It remains to note that the Schur 
complement S in (49) is symmetric positive definite because S = 5rr — □ 

Theorem 4 allows us to use the conjugate gradient method for the iterative solution of (49). In the 
next section, we describe the BDDC method used as a preconditioner for S. 


5. THE BDDC PRECONDITIONER 

In this section, we formulate the BDDC method for the solution of (49). The algorithm can be 
viewed as a generalisation of [6]. However, we follow the original description from [11], which 
better reflects our implementation. 

One step of BDDC provides a two-level preconditioner for the conjugate gradient method applied 
to solving problem (49). It is characterised by the selection of certain coarse degrees of freedom 
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based on primary degrees of freedom at interface F. The main coarse degrees of freedom in this 
paper are arithmetic averages over faces, defined as subsets of degrees of freedom shared by the 
same two substructures. In addition, corner coarse degrees of freedom, defined as any selected 
Lagrange multiplier at the interface, are used. Substructure edges, defined as subsefs of degrees of 
freedom shared by several substructures, may also appear (see Remark 5). 

The BDDC method introduces constraints which enforce continuity of functions from Ar at 
coarse degrees of freedom among substructures. This gives raise to the space Ar, which is given as 
the subspace of Ar of functions which satisfy these continuity constraints. In particular, 

Ar C Ar C Ap. (52) 


Remark 5 

In three spatial dimensions, several triangular elements can be connected at a single Lagrange 
multiplier in a star-like conhguration (cf. Fig. 1). A similar statement holds for line elements 
considered in two- or three-dimensional space. This fact may lead to the presence of substructure 
edges and even vertices (defined as degenerate edges consisting of a single degree of freedom), and 
we may prescribe also edge averages as constraints. As mentioned above, we also select corners 
as coarse degrees of freedom. While essentially any degree of freedom at the interface T can be a 
comer, we select them by the face-based algorithm [21]. This algorithm considers all vertices as 
corners, and, in addition, it selects three geometrically well distributed degrees of freedom from the 
interface between two substructures sharing a face into the set of comers. Although considering 
corners is not the standard practice with RTq finite elements, in our experience, comers improve 
convergence for numerically difficult problems, as can be observed for the engineering applications 
presented in Section 8. 

We now proceed to the formulation of operators used in the BDDC method. The choice of 
constraints determines the construction of matrices DU Each row of D* defines one coarse degree of 
freedom at substructure D®, e.g., a comer corresponds to a single 1 entry at a row and an arithmetic 
average to several I’s at a row. The coarse basis functions one per each substructure coarse 
degree of freedom, are computed by augmenting the matrices from (43) with D®, and solving the 
augmented systems with multiple right-hand sides 


■ A® 

rzT 

TjiT 

Bf,i 


0 ■ 


- - 


■ 0 - 

B® 

-C® 

CiT 

CiT 

0 


Z® 


0 


-c>./ 

-C}i 

CliT 

Lyj 

0 



= 

0 


-Cfr 

—C^i 

Orr 

^®T 




0 

0 

0 

0 

D® 

0 


A* 


. I _ 


1 = 1,...,Ns, (53) 


where I is the identity matrix, and AT®, Z®, and are auxiliary matrices not used any further. As 
shown in [35], the local coarse matrix is obtained as a side product of solving (53) as 


Sfc = [A®^ Z®^ 4>f 


A® 

B^t 

A® 

-C® 

Bfj 

-Cf,i 

Bf.r 



jDiT 

Bf,i 

r}iT 

B^.r 

/^iT 

C^iT 

-cji 

CiT 

—fri 

—Cyj 

Orr 



- - 


Z® 



- 1 

. . 


-L\ 


(54) 


Let us define, similarly to (47), operators Rf. that relate vectors of local coarse degrees of 
freedom to the vector of global coarse degrees of freedom fj,c as 


dC — RcdC- 


(55) 


The global coarse matrix Sec is then obtained by assembling the local contributions as 


Ns 

See = Y.^c Shelve- 


2=1 


(56) 
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Finally, let us define the scaling operators 

i = (57) 

which are given as diagonal matrices of weights that satisfy 


Ns 

2 = I. (58) 

Z=1 

More details on the selection of diagonal entries in are given in Section 6. 

With this selection of spaces and operators, we are ready to formulate the BDDC preconditioner. 

Algorithm 6 

The BDDC preconditioner Mbddc : rr S Ar —Ar € Ar is defined in the following steps; 


1. Compute the local residuals 


4 = VF*i?Vr, 


(59) 


2. Compute the substructure corrections by solving the local Neumann problems 


■ A* 

B^t 

tdiT 

RtT 

0 ■ 


at* 


■ 0 - 

S* 

-C* 

r^iT 

riiT 

0 


z* 


0 

Byi 


W/ 

r<iT 

Lyj 

0 



= 

0 

Br,r 


Op/ 

Opr 

DiT 




4 

0 

0 

0 

D* 

0 


^* 


0 


3. Compute the coarse correction by collecting the coarse residual 


(60) 


Ns 

solving the global coarse problem 

Sec VC =rc, 

and distributing the coarse correction 

p^rc = ^rRhvc. i = l,....Ns. 
4. Combine and average the corrections 


Ar 


Ns 


2 = 1 


+ Vre)■ 


( 61 ) 


(62) 

(63) 


(64) 


We note that the factorisations of the matrices from (53) are also used for each solution of (60). 
In order to apply the existing BDDC theory for elliptic problems (e.g. [13, 36]) to the proposed 
preconditioner, we introduce some additional notation and make a few observations. Due to (44), 
the substructure corrections in (60) can be written equivalently as 


-S'* 



Vvis 




r 


0 


and the construction of coarse basis functions in (53) as 


-S* 



■ $*, ■ 


■ 0 ■ 


L* 


I 


(65) 


(66) 
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Next, let us formally write the operators and vectors in the block form 



■ Ai ■ 


7?i 


■ 


■ si 

Ar = 


,77 = 


,w = 


,S = 



_ _ 


rNs 




SNs 


. (67) 


By grouping the steps of Algorithm 6 and using (65), the operator of the BDDC preconditioner 
can be formally written as 

Mbddc = R^WS-^WR, (68) 



diag 


[4^ 


_S^ J^iT 
0 


-1 


0 


'Ns 


■Ns 


Vi=l / Vi=l 

(69) 


iT 

r 


The first term in S~^ corresponds to substructure corrections and the second term to the coarse 
correction (steps 2 and 3 of Algorithm 6), and If^i is the identity matrix in Ap. From (68) and (69), 
one can readily see that Mbddc is symmetric. 


Assumption 7 
Let us assume that 


nulls'* _L nullD*. 


(70) 


In order to satisfy Assumption 7, we must prescribe enough coarse degrees of freedom as 
constraints along with Robin boundary conditions (19)-(20) and (22) at each fracture within 
substructure fl*. Since constrains in D* are linearly independent, has full column rank. In 
particular. Assumption 7 is satisfied when arithmetic averages are used on each substructure face 
(and eventually edge) as constraints. 

Lemma 8 

The operator S~^ in preconditioner Mbddc is symmetric and positive definite on the space Ap. 


Proof 


The space Ar is decomposed into the substructure spaces and the coarse space, 

Ar = ArA © Arc- (71) 

To achieve this splitting, each local space Ap is decomposed into subspaces 

Nf = null D* 0 range (72) 

corresponding to the substructure space Ap^, and the coarse space on substructure fl*. Arcing 
respectively. To analyse this decomposition, let us recall that S'* is a positive semi-definite matrix 
and write (66) in detail as 


79*$^ = I. 


From (73), 


range (S*T>r) C range 79*^ _L null79\ 

which in turn, similarly to [36, Lemma 8], gives for any xa ^ ?c ^ range $p 


X^lS% = 0 . 


(73) 

(74) 

(75) 

(76) 
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From (74), the matrix $p has full column rank, and 

null 77* n range ‘hp = 0. (77) 


Finally, from Assumption 7 and (72), 

null S’' C range 4>p. 


(78) 


Decomposition of the subdomain space (72) implies decomposition of a function G Ap to 
C* = Ca + Cc’ where G null 77®, e range $p, and (’a = 0 by (76). 

Let us first analyse the substructure corrections. Following [3, Section 3.3], the matrix of (65) is 
invertible due to Assumption 7. If we define, in addition, a matrix Q’ with orthonormal columns 
forming a basis of null 77®, i.e. 

range Q® = null 77®, = I, (79) 


we have 


[ 


hi 0 


The matrix S’Q’) ^ is symmetric positive definite, and consequently, for any (’ G Ap, 


■ -S’ 

D’T ■ 

-1 

f 1 

77® 

0 


0 


= -Q® {Q’^S’Q’) Q’^. 


(80) 


- {Q’^S’Q’) ^ Q’^C < 0 (81) 

with equality if g range <i)p. 

Next, let us turn towards the coarse correction. Formula (54) for S’(j(j can be written equivalently 
as 

S’cc = -^f S’iS>’r- ( 82 ) 

Since the term on the right-hand side is just a (negative) Galerkin projection of the positive semi- 
definite matrix S’, matrix is symmetric negative semi-definite. If at least one substructure is 
equipped with natural boundary conditions, the matrix Sec assembled by (56) becomes symmetric 
negative dehnite, and so is Sqq. 

We have just verified the negative definiteness of the principal parts of S~^, and the desired 
positive definiteness is obtained through the change of sign in front of the braces in (69). □ 


In view of Lemma 8, the standard condition number bound follows from [13, Lemma 2]. 
Theorem 9 

The condition number k of the preconditioned operator MbddcS satisfies 

||7?7?'^lLAr 

K < u) = max -2— 

ArGAr ll-^rUg 

The norms in (83) are induced by the matrix S defined in (67) for all functions Ar G Ap. 

In addition, in the case of a single mesh dimension in either 2D or 3D, and under the assumption 
of substructure-wise constant hydraulic conductivities, it has been also derived in [6, Lemma 5.5 
and Theorem 6.1] that the condition number bound uj satisfies 

t^<c(^l + log|^^ , (84) 

where H is the characteristic size of geometric substructures. We note that the bound (84) implies 
that, for a fixed relative subdomain size H/h, the condition number is independent of the problem 
size. 

It is worth emphasising that Theorem 9 is also valid for combined mesh dimensions. Flowever, 
several simpliheations are employed in [6] to obtain (84), which are not satisfied in the set-up 
considered in this paper. In particular (i) hydraulic conductivity coefficient here is, in general, not 
substructure-wise constant nor isotropic, (ii) it is not clear whether, in presence of fractures, the 
interpolation operator onto a conforming mesh introduced in [6] can be constructed and bounded in 
the norm. 


2 

(83) 
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6. SCALING WEIGHTS IN BDDC 


Let us now discuss the choice of entries in the diagonal weight matrices IL*. These matrices play 
an important role in the BDDC method, both in the theory, cf. Theorem 9 or [6, 13, 37], and in 
the computations, cf. [19, 38]. Three possible choices are also studied numerically in Section 8.2. 
The basic choice is presented by the arithmetic average taken from values at the neighbouring 
substructures. In this simplest construction, the entry corresponding to Lagrange multiplier Ap ^ is 
given by the inverse counting function as 


W = - 

cardiXj 


(85) 


where card{Xj) is the number of substructures in the setlj of indices of substructures to which Ap ^ 
belongs. Lor 2D or 3D meshes without fractures, VLU = 1/2 for Raviart-Thomas elements. However, 
since several two-dimensional fractures can meet in our setting, smaller weights can occasionally 
appear at such regions. 

While arithmetic average is sufficient for problems with homogeneous coefficients, it is well 
known that for problems with large variations in material properties along the interface, it is 
necessary to incorporate their values into the (weighted) average to obtain a robust method. This 
gives rise to the p-scaling, for which 


Wl 


Pi 

E cardlXi) 

k=l Pk 


( 86 ) 


where pk is a material characteristic for substructure This choice is robust with respect to jumps 
in coefficients across the interface, cf. [6, 9]; however, coefficients are assumed constant for each 
substructure. This requirement is very restrictive for practical computations with quickly varying 
coefficients, and we employ a generalisation which takes into account the material coefficient of the 
element to which the Lagrange multiplier Ap^ corresponds. In our case, we use pi = d/tr(Ik“^), 
where d G {1,2,3} is the dimension of the element T*. This value can be seen as a representative 
hydraulic conductivity on the element. 

Linally, we propose a modification of the popular scaling by diagonal stiffness [19]. In the 
usual diagonal stiffness approach, the optimal weight, which is the diagonal entry of the Schur 
complement, is approximated by the diagonal entry of the original substructure matrix. However, 
this is not directly applicable to the indefinite system (43), as, in general, matrix C® contains 
only seldom nonzeros on the diagonal. Lor this reason, we approximate the diagonal of the Schur 
complement as 


Wh=ChT,n 


Ai 


kk 


(87) 


where the index k corresponds to the row in block A® of the element face to which the Lagrange 
multiplier Ap ^ belongs. 

Using the diagonal stiffness scaling in connection with the standard Lagrange finite elements may 
lead to poor convergence for problems with rough interface [19, 38], for which the diagonal stiffness 
can vary quickly even for smooth problems with constant coefficients on uniform meshes. This is a 
severe issue for practical computations, in which graph partitioners are typically used for creating 
substructures. However, this issue is not as pronounced for Raviart-Thomas elements, for which 
only one element contributes to the stiffness on the diagonal at an interface degree of freedom, and 
thus irregularities caused by changing number of elements contributing to an interface weight cannot 
occur. On the other hand, an advantage of the diagonal stiffness scaling is the fact, that—unlike the 
p-scaling—it takes into account the shape and relative sizes of elements, which vary considerably 
in engineering applications, as well as the effect of 5d introduced in (29) and (36). Unless stated 
otherwise, scaling (87) is used in the computations presented in Section 8. 
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7. THE PARALLEL SOLVER 

The basis for an efficient parallel implementation of the method described in previous sections was 
obtained by combining two existing open-source software packages: the finite element package 
Flowl23d^ (version 1.6.5) for underground fluid flow simulations and the BDDC-based solver 
BDDCMLf (version 2.0) used for the solution of the resulting system of equations. However, minor 
changes have been made to both codes to support the specific features, such as the weights (86) 
and (87). 

The Flowl23d package has been developed for modelling complex behaviour of underground 
water flow and pollution transport. However, only the simple flow in a fully saturated porous 
media described by Darcy’s law is considered in this paper. To accurately account for fractures 
in the medium, such as granite rock, the solver allows us to combine finite elements of different 
dimensions: the three-dimensional elements of porous media are combined with two-dimensional 
elements modelling planar fractures, which may be in turn connected in one-dimensional elements 
for channels. Raviart-Thomas elements are consistently used throughout such discretisation. 
Although the fractures are also modelled as porous media, their hydraulic conductivity is by 
orders of magnitude higher than that of the main porous material of the domain. In addition, the 
finite element discretisations are typically not uniform within the domain, and the relative sizes 
of elements may also vary by orders of magnitude. Both these aspects give rise to very poorly 
conditioned linear systems, which are very challenging for iterative solvers. The Flow 123d solver 
has been developed for over 10 years and it is written in C/C++ programming language with object- 
oriented design and parallelism through MPI. 

The BDDCML is a library for solving algebraic systems of linear equations by means of the 
BDDC method. The package supports the Adaptive-Multilevel BDDC method [22] suitable for very 
high number of substructures and computer cores, although we only use the standard (non-adaptive 
two-level) BDDC method from [6, 11] for the purpose of this paper. The BDDCML library is 
typically interfaced by finite element packages, which may provide the division into substructures. 
This feature is used in our current implementation, in which the division into non-overlapping 
substructures is constructed within the Flowl23d using the METIS (version 5.0) package [39]. 
One substructure is assigned to a processor core in the current set-up of the parallel solver, although 
BDDCML is more flexible in this respect. The library performs the selection of additional corners 
by the face-based algorithm from [21]. The BDDCML package is written in Lortran 95 and 
parallelised through MPI. 

The BDDCML solver relies on a serial instance of the MUMPS direct solver [40] for the solution 
of each local discrete Dirichlet problem (45) as well as for the solution of each local discrete 
Neumann problem (60). The coarse problem (62) is solved by a parallel instance of MUMPS. The 
main difference from using BDDCML for symmetric positive definite problems is the need to use 
the LDL"’^ factorisation of general symmetric matrices for problems (45), which are saddle-point 
(i.e. indefinite) systems in the present setting. 

Although the original system (37) is indefinite, system (49) is symmetric positive definite, which 
allows the use of the preconditioned conjugate gradient (PCG) method. One step of BDDC is used 
as the preconditioner within the PCG method applied to problem (49). The matrix of problem (49) 
is not explicitly constructed in the solver, and only its actions on vectors are computed following 
(45)-(48). 

Remark 10 

In our implementation, we change the sign neither in the action of S'* (46) nor in the action of the 
preconditioner Mbddc (64). Since both are then strictly negative definite, the product MbddcS is 
the same as if both signs were changed, and the PCG method runs correctly. In this way, no changes 
are necessary in an implementation developed for symmetric positive definite problems. 


thttp://flowl23d.github.io 

thttp://users.math.cas.cz/~sistek/software/bddcmi.htmi 
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8. NUMERICAL RESULTS 

In this section, we investigate the performance of the algorithm and its parallel implementation 
on several benchmark problems in 2D and 3D, and on two geoengineering problems of existing 
localities in 3D. Eor the two benchmark problems without fractures we perform weak scaling tests. 
Eor the benchmark problem with fractures and for the geoengineering problems, we perform strong 
scaling tests with the problem size fixed and increasing number of processor cores. In all cases, the 
PCG iterations are run until the relative norm of residual ||/ll^|| < 10“^. If not stated otherwise, 
the proposed scaling by diagonal stiffness (87) is used within the averaging operator of BDDC. 


8.1. Results for benchmark problems 

Eirst, the performance of the solver is investigated on a unit square and a unit cube discretised 
solely using two-dimensional and three-dimensional finite elements, respectively. Eor this reason, 
block C in system (37), which is related to combining elements of different dimension, is zero, 
and the problem reduces to the standard problem (14). The sequence of unstructured meshes is 
approximately uniform for both problems, and the problems do not contain any jumps in material 
coefficients. In Eigs. 2 and 3, example meshes and the resulting pressure head and velocity fields 
are presented. While gravity is present in the 3D case, its effect is not considered in the 2D case. 

The results of the weak scaling tests are summarised in Tables I and II. To give a better view, the 
resulting solution times for different problem sizes are also visualised in Eig. 4. In these tables, N 
denotes the number of substructures and processors, n is the size of the global problem (37), up is 
the size of the interface problem (49), n/ denotes number of faces, Uc denotes number of comers, 
‘its.’ stands for resulting number of PCG iterations, and ‘cond.’ is the approximate condition number 
computed from the Lanczos sequence in PCG. We report separately the time spent in preconditioner 
set-up, the time spent by PCG iterations, and the total time for the whole solve. 

In these weak scaling tests, the number of unknowns per core is kept approximately constant 
around 10®. These weak scaling tests were performed using up to 64 cores of the SGI Altix UV 
supercomputer at the Supercomputing Centre of the Czech Technical University in Prague. The 
computer contains twelve Intel Xeon processors, each with six cores at frequency 2.67 GHz. Intel 
compilers version 12.0 were used. 

The numbers of PCG iterations and condition number estimates in Tables I and II confirm the 
expected numerical scalability of the BDDC method, which is well known for symmetric positive 
definite problems as well as for Darcy’s flow problems [6, 37]. The slight irregularities in the 
condition number in Table II are probably caused by using non-nested unstructured meshes. 

Looking at times in these tables and in Eig. 4, we can see almost optimal scaling, with only mild 
growth of times with number of cores. The numbers of PCG iterations are higher in 3D, and the 
time spent in PCG iterations grows proportionally, while the time spent in the set-up phase does 
not differ considerably between two-dimensional and three-dimensional setting and dominates the 
overall time. 


N 

n 

n/N 

nr 


Tic 

its. 

cond. 

time (sec) 

set-up PCG solve 

2 

207k 

103k 

155 

1 

2 

7 

1.37 

8.3 

1.6 

9.9 

4 

440k 

110k 

491 

5 

10 

8 

1.60 

12.2 

2.2 

14.4 

8 

822k 

103k 

1.2k 

13 

26 

9 

1.78 

11.0 

2.5 

13.5 

16 

1.8M 

111k 

2.8k 

33 

66 

8 

1.79 

14.3 

2.7 

17.0 

32 

3.3M 

104k 

5.9k 

74 

148 

9 

1.79 

12.1 

3.3 

15.4 

64 

7.2M 

113k 

13.0k 

166 

332 

9 

1.85 

14.8 

4.4 

19.2 


Table I. Weak scaling test for the 2D square problem, each substructure problem contains approx. 100k 

unknowns. 
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pressure 

0.974512 

10.8 

iO.4 

*0 

-0,4 

t-0,8 

-0,974142 



velocity mag. 

13.12612 

■ 

f 

2.629e-5 


Figure 2. Example of solution to the model square problem containing only 2D elements, plot of pressure 

head with mesh (left) and velocity vectors (right). 



pressure 

1.821991 


r' 

-1.857597 



Figure 3. Example of solution to the model cube problem containing only 3D elements, plot of pressure 

head with mesh (left) and velocity vectors (right). 


N 

n 

n/N 

nr 


Tie 

its. 

cond. 

time (sec) 

set-up PCG solve 

2 

217k 

108k 

884 

1 

3 

11 

2.88 

11.7 

2.3 

14.0 

4 

437k 

109k 

2.3k 

6 

18 

12 

3.04 

11.7 

2.5 

14.2 

8 

945k 

118k 

5.7k 

21 

63 

15 

12.00 

15.4 

4.0 

19.3 

16 

1.6M 

103k 

12.8k 

56 

168 

16 

6.58 

12.9 

4.0 

17.0 

32 

3.4M 

106k 

29.8k 

132 

401 

18 

10.10 

15.4 

5.2 

20.6 

64 

6.1M 

95k 

59.6k 

307 

931 

19 

16.58 

13.7 

6.3 

20.0 


Table II. Weak scaling test for the 3D cube problem, each substructure problem contains approx. 100k 

unknowns. 


The next benchmark problem is considerably more complicated. It consists again of a unit cube, 
which now contains four planar fractures aligned with diagonals of a 2D cross-section. These four 
planar fractures meet at a ID channel in the centre of the cross-section. Therefore, the problem 
contains the full possible combination of 3D, 2D and ID finite elements. The tensor k is isotropic, 
thus it is just a scalar multiple of identity. The corresponding scalar value is set to 10, 1, and 0.1 for 
ID, 2D, and 3D elements, respectively. 

We perform a strong scaling test with this problem, keeping the mesh size fixed with 
approximately 2.1 million elements and 14.6 million degrees of freedom. In Fig. 5, the 
computational mesh and the resulting pressure head and velocity fields are presented. This scaling 
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a 



Figure 4. Weak scaling test for the 2D square problem (left), and the 3D cube problem (right), approx. 100k 
unknowns per core. Computational times separately for set-up and PCG phases, and their sum (solve). 


test was computed on the Cray XE6 supercomputer Hector at the Edinburgh Parallel Computing 
Centre. This supercomputer is composed of 2816 nodes, each containing two AMD Opteron 
Interlagos processors with 16 cores at 2.3 GHz. GNU compilers version 4.6 were employed. 

Results of the strong scaling test are summarised in Table 111, and the computing times are 
visualised also in Eig. 6 together with the parallel speed-up. The reference value for computing 
speed-up is the time on 16 cores, and the speed-up on np processors is computed as 


^np — 


16 fi6 

^np 


( 88 ) 


where is the time on np processors. 

We can see that the number of PCG iterations grows with the number of substructures for this 
problem which is also confirmed by the growing condition number estimate. While the time spent in 
set-up phase scales very well, the time spent in PCG grows together with the number of iterations. 
The reason for this growth seem to be related to the larger interface, at which more numerical 
difficulties appear. This seems to be related to more 1D-2D and 2D-3D connections at the interface 
and makes this difficult problem a good candidate for using the Adaptive BDDC method [22, 41]. 
However, this will be the subject of a separate study. 



pressure 

0.986152 

P0,8 

U.4 

to 

i-0.4 

1 - 0.8 

-0,985993 



velocity mag. 
0.369957 

*0,3 
^ 0,2 
jO.l 

0.02923 


Figure 5. Example of solution to the model cube problem containing ID, 2D, and 3D elements, plot of 
pressure head with mesh and fractures (left) and velocity vectors (right). 
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N 

n/N 

nr 


Tic 

its. 

cond. 

time (sec) 

set-up PCG solve 

16 

912k 

47k 

53 

159 

26 

59.3 

171.6 

84.5 

256.2 

32 

456k 

65k 

126 

380 

48 

2091.0 

90.1 

109.8 

200.0 

64 

228k 

86 k 

301 

914 

81 

1436.1 

36.8 

77.1 

114.0 

128 

114k 

116k 

689 

2076 

109 

2635.8 

14.3 

43.1 

57.4 

256 

57k 

151k 

1436 

4365 

164 

1700.5 

6.7 

31.2 

38.0 

512 

28k 

196k 

3021 

9244 

254 

42614.5 

4.0 

26.9 

30.9 


Table III. Strong scaling test for the cube problem with ID, 2D, and 3D elements, size of the global problem 

is n = 14.6M unknowns. 



number of processors 



Figure 6. Strong scaling test for the cube problem with ID, 2D, and 3D elements and 14.6M unknowns, 
computational time (left) and speed-up (right) separately for set-up and PCG phases, and their sum (solve). 


8.2. Results for geoengineering problems 

The performance of the algorithm and its parallel implementation has been investigated on two 
engineering problems of underground flows within real geologic locations. For both problems, the 
porous medium is fractured granite rock, with the fractures modelled by two-dimensional elements. 

The first problem is the Melechov locality, which models one of the candidate sites for a 
nuclear waste deposit to be build within the Czech Republic in future. The goal is to model the 
underground flow and estimate the speed at which an eventual radioactive pollution would spread. 
The computational mesh contains 2.1 million finite elements resulting in 15 million unknowns. The 
geometry of the problem with the resulting distribution of piezometric head and the finite element 
mesh is presented in Fig. 7. The problem contains vertical two-dimensional fractures visualised in 
Fig. 8 . The maximal hydraulic conductivity within the fractures is 6.3-10^ ms“^, while the minimal 
conductivity of the outer material is 6 . 0 - 10 “^ ms“^, the transition coefficient 0 - 3=1 s“^, and the 
effective thickness of fractures ^2 = 0.1 m. 

We perform a strong scaling test for this problem, keeping the problem size fixed and increasing 
the number of substructures and computing cores. An example of division into 64 substructures is 
presented in Fig. 8 . The scaling test was computed on the Hector supercomputer. 

Table IV summarises the results of this test. We can still see some growth of the number of 
iterations with the number of substructures, which is however much milder than the growth observed 
for the unit cube with fractures in Table III. Correspondingly, the times reported in Table IV and 
visualised in Fig. 9 show an optimal scaling of the solver over a large range of core counts. 

The second engineering model is the locality around the Bedfichov tunnel. The main purpose of 
this 2.1 km long tunnel near the city of Liberec in the north of the Czech Republic to accommodate 
water pipes, which supply the city by drinking water from a reservoir in the mountains. Flowever, 
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Figure 7. The problem of the Melechov locality containing 2D and 3D elements, mesh contains 2.1M 
elements and 15M unknowns. Plot of the piezometric head. Data by courtesy of Jifina Kralovcova. 



Figure 8. The problem of the Melechov locality, the system of fractures (left) and an example division into 

64 substructures (right). 


N 

n/N 

Uy 


Tie 

its. 

cond. 

set-up 

time (sec) 
PCG 

solve 

16 

934k 

36k 

32 

96 

40 

53.0 

131.4 

144.1 

275.6 

32 

467k 

54k 

76 

228 

70 

878.3 

47.5 

112.9 

160.4 

64 

233k 

82k 

186 

561 

67 

202.4 

17.4 

50.2 

67.7 

128 

117k 

116k 

528 

1592 

69 

237.6 

7.9 

23.1 

31.1 

256 

58k 

155k 

1235 

3747 

96 

5577.0 

4.0 

14.7 

18.8 

512 

29k 

207k 

2699 

8256 

106 

1658.1 

2.2 

8.3 

10.5 

1024 

15k 

271k 

5711 

17581 

119 

11554.5 

2.1 

7.0 

9.2 


Table IV. Strong scaling test for the problem of the Melechov locality containing 2D and 3D elements, size 

of the global problem is n = 15M unknowns. 


this locality is also a valuable site for experimental geological measurements performed inside the 
tunnel. 

The model aims at describing the flow in the granite rock surrounding the tunnel. The 
computational mesh consists of 1.1 million elements leading to 7.8 million unknowns. The mesh 
with the plot of resulting piezometric head is presented in Fig. 10. The system of fractures and an 
example division into 256 substructures are visualised in Fig. 11. The hydraulic conductivity of the 
fractures is 10 “^ ms“^, while the conductivity of the outer material is 10 “^° ms“^, the transition 
coefficient 0-3 = 1 s“^, and the effective thickness of fractures ^2 = 1-1 m. 
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number of processors 



Figure 9. Strong scaling test for the problem of the Melechov locality containing 2D and 3D elements and 
15M unknowns, computational time (left) and speed-up (right) separately for set-up and PCG phases, and 

their sum (solve). 


Although the mesh contains fewer finite elements than the one of the Melechov locality model, 
this problem is considerably more complicated. This is caused mainly by the presence of relatively 
very small and irregularly shaped finite elements in the vicinity of the tunnel and near the cross- 
sections of fractures (see Fig. 12) generated by the mesh generator. 

The results of a strong scaling test are summarised in Table V. As before, the times are also plotted 
in Fig. 13. Although the number of iterations is not independent of the number of substructures, the 
growth is still small. Consequently, the computing times, and especially the time for set-up, scale 
very well over a large range of numbers of substructures. The observed super-optimal scaling may 
be related to faster factorisation of the smaller local problems by the direct solver for indefinite 
matrices. 

Table VI summarises an experiment performed to analyse the effect of using comers in the 
construction of the coarse space in BDDC. As has been mentioned in Section 5, using Raviart- 
Thomas finite elements does not lead to ‘natural’ comers as cross-points shared by several 
substructures. On the other hand, the notion of comers was generalised to any selected interface 
degree of freedom, at which continuity of functions from the coarse space is required. Such 
generalisation is important for the well-posedness of the local problems for unstmctured meshes e.g. 
in elasticity analysis [21]. This is also the default option for BDDCML, in which selection of corners 
is performed at each face between two substructures. Adding comers improves the approximation 
properties of the coarse space at the cost of increasing the size of the coarse problem. Table VI 
compares the convergence for variable number of substructures without and with constraints at 
corners. Column ‘with corners’ in Table VI corresponds to results in Table V, and it is repeated for 
comparison. We can see that while the effect of corners on convergence is small for smaller number 
of substructures, the improvement of the coarse problem and the approximation of BDDC becomes 
more significant for higher numbers of cores. Looking at times in Table VI, the additional time spent 
in the set-up phase due to higher number of constraints when using comers is compensated by the 
lower number of PCG iterations, resulting in lower overall times. Thus, using additionally selected 
comers appears beneficial for complicated engineering problems like this one. 

In the final experiment, we compare the effect of different averaging techniques on the 
convergence of BDDC. In Table VII, results of the strong scaling test for arithmetic averaging (85), 
the modified p-scaling (86), and the proposed scaling by diagonal stiffness (87) are summarised. 
The final column corresponds to the results from Table V, which are repeated here for comparison. 

Table VII suggests, that while the simple arithmetic averaging does not lead to satisfactory 
convergence for this problem, the modified p-scaling and the diagonal scaling mostly lead to similar 
convergence. However, while the former provides slightly better convergence for several cases, it 
also leads to irregularities for certain divisions, for which the BDDC method with this averaging 
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Figure 10. The Bedfichov tunnel problem containing 2D and 3D elements, mesh contains I.IM elements 
and 7.8M unknowns. Plot of the piezometric head. Data by courtesy of Dalibor Frydrych. 



Figure 11. The Bedfichov tunnel problem; the system of fractures (left) and an example division into 256 

substructures (right). 



Figure 12. Example of difficulties in the mesh of the Bedfichov tunnel problem; detail of the tunnel geometry 
with fine elements (left) and enforced refinement at an intersection of fractures (right). 


converges rather poorly. Therefore, the proposed scaling (87) can be recommended as the most 
robust choice among the three tested options. 
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N 

n/N 

nr 

«/ 

Tie 

its. 

cond. 

set-up 

time (sec) 
PCG 

solve 

32 

245k 

20k 

106 

322 

112 

1514.1 

110.3 

144.0 

254.3 

64 

123k 

28k 

192 

597 

63 

117.7 

42.2 

36.0 

78.3 

128 

61k 

45k 

413 

1293 

75 

194.4 

13.4 

16.8 

30.3 

256 

31k 

72k 

902 

2791 

119 

526.7 

4.2 

10.9 

15.1 

512 

15k 

110k 

2009 

6347 

137 

1143.4 

1.8 

7.1 

9.0 

1024 

8k 

155k 

4575 

14725 

173 

897.0 

1.6 

8.0 

9.7 


Table V. Strong scaling test for the problem of the Bedfichov tunnel containing 2D and 3D elements, size 

of the global problem is n = 7.8M unknowns. 




number of processors number of processors 


Figure 13. Strong scaling test for the problem of the Bedfichov tunnel containing 2D and 3D elements and 
7.8M unknowns, computational time (left) and speed-up (right) separately for set-up and PCG phases, and 

their sum (solve). 


N 

n/N 

its. 

without corners 

, time (sec) 

cond. ^ 

set-up PCG 

solve 

its. 

cond. 

with comers 

time (sec) 
set-up PCG 

solve 

32 

245k 

131 

1789.6 

107.5 

175.0 

282.5 

112 

1514.1 

110.3 

144.0 

254.3 

64 

123k 

70 

122.2 

40.3 

40.4 

80.7 

63 

117.7 

42.2 

36.0 

78.3 

128 

61k 

96 

208.5 

10.9 

21.6 

32.6 

75 

194.4 

13.4 

16.8 

30.3 

256 

31k 

139 

541.9 

3.7 

12.5 

16.2 

119 

526.7 

4.2 

10.9 

15.1 

512 

15k 

197 

1418.9 

1.4 

10.0 

11.4 

137 

1143.4 

1.8 

7.1 

9.0 

1024 

8k 

312 

3779.4 

1.0 

14.5 

15.6 

173 

897.0 

1.6 

8.0 

9.7 


Table VI. Effect of using corners. Problem of the Bedfichov tunnel containing 2D and 3D elements, size of 
the global problem is n = 7.8M unknowns. In column ‘without corners’, no additional corners are selected 
in BDDC. In column ‘with corners’, additional comers are selected. 


9. CONCLUSION 

A parallel solver for the mixed-hybrid finite element formulation based on Darcy’s law has been 
presented. The software combines an existing package Flowl23d developed for problems in 
geophysics with BDDCML, a parallel library for solution of systems of linear algebraic equations 
by the BDDC method. 

In geoengineering applications, the mathematical model is applied to geometries with the 
presence of fractures. In the present approach, the flow in these fractures is also modelled by Darcy’s 
law, although the hydraulic conductivity of the porous media is considered by orders of magnitude 
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N 

n/N 

nr 

nf 

ric 

arithmetic avg. 
its. cond. 

mod 

its. 

. p-scal. 
cond. 

diagonal seal, 
its. cond. 

32 

245k 

20k 

106 

322 

637 

9811.7 

no 

1467.8 

112 

1514.1 

64 

123k 

28k 

192 

597 

618 

10254.1 

62 

115.1 

63 

117.7 

128 

61k 

45k 

413 

1293 

2834 

l.Oe-i-11 

206 

401641.4 

75 

194.4 

256 

31k 

72k 

902 

2791 

799 

11172.9 

117 

512.9 

119 

526.7 

512 

15k 

110k 

2009 

6347 

883 

15449.6 

136 

1160.1 

137 

1143.4 

1024 

8k 

155k 

4575 

14725 

n/a 

2.5e-i-10 

504 

99023.6 

173 

897.0 


Table VII. Comparison of different averaging techniques for the Bedfichov tunnel containing 2D and 3D 
elements, size of the global problem is n = 7.8M unknowns. 


higher. These fractures are modelled by finite elements of a lower dimension. In the discretised 
model, ID, 2D and 3D finite elements are coupled together through Robin boundary conditions. 
These coupling terms lead to a modification of the usual saddle-point matrix of the system, in which 
a new non-zero block appears on the diagonal. 

The BDDC method is employed for the solution of the resulting system of linear algebraic 
equations. BDDC is based on iterative substructuring, in which the problem is first reduced to 
the interface among substructures. The Schur complement is not built explicitly. Instead, only 
multiplications by the matrix are performed through solving a discrete Dirichlet problem on each 
substructure. In the setting of the mixed-hybrid problem, the interface is built only as a subset of 
the block of Lagrange multipliers, while remaining unknowns belong to interiors of substructures. 
Although the original problem is symmetric indefinite, the system reduced to the interface is 
symmetric positive definite. This is also shown to hold for the case with fractures in the present 
paper. Consequently, the PCG method is used for the solution of the reduced problem. However, 
unlike the symmetric positive definite problems, a direct solver for symmetric indefinite matrices 
needs to be used for the factorisation and repeated solution of local problems on substructures. 

One step of the BDDC method is used as the preconditioner for the PCG method run on the 
interface problem. A modification of the diagonal stiffness scaling has been introduced. It is 
motivated by difficult engineering problems, for which it performs significantly better than other 
two applicable choices—the arithmetic averaging and the modified p-scaling. Arithmetic averages 
over faces between substructures are used as the basic constraints defining the coarse space. In 
addition, corners are selected from unknowns at the interface using the face-based algorithm. While 
comers are not required by the theory, they are shown to improve both the convergence and the 
computational times for complicated problems. 

The performance of the resulting solver has been investigated on three benchmark problems in 
2D and 3D. Both weak and strong scaling tests have been performed. On benchmark problems with 
single mesh dimension, the expected optimal convergence independent of number of substructures 
has been achieved. Correspondingly, the resulting parallel scalability has been nearly optimal for 
the weak scaling tests up to 64 computer cores. 

The strong scaling tests were presented for a benchmark problem of a unit cube and for two 
engineering problems containing large variations in element sizes and hydraulic conductivities, 
using up to 1024 computer cores and containing up to 15 million degrees of freedom. The 
convergence for the unit cube problem with all three possible dimensions of finite elements slightly 
deteriorated by using more substructures, and this translated to sub-optimal parallel performance. 
However, for the two engineering applications, in which only 3D and 2D elements are combined, 
the BDDC method has also maintained good convergence properties with the growing number of 
substructures, resulting in optimal, or even super-optimal parallel scalability of the solver. It has been 
also shown that the proposed modification of the diagonal stiffness scaling plays an important role 
in achieving such independence for the challenging engineering problems presented in the paper. 


ACKNOWLEDGEMENT 


Preprint version. 




26 


J. SISTEK, J. BREZINA, B. SOUSEDIK 


The authors are grateful to Jifina Kralovcova and Dalibor Frydrych for providing the geoengineering models. 


REFERENCES 

1. Gulbransen AF, Hauge VL, Lie KA. A multiscale mixed finite-element method for vuggy and naturally-fractured 
reservoirs. SPE Journal 2010; 15(2):395^03. 

2. Martin V, Jaffre J, Roberts JE. Modeling fractures and barriers as interfaces for flow in porous media. SIAM Journal 
on Scientific Computing 2005; 26(5):1667-1691. 

3. Benzi M, Golub GH, Liesen J. Numerical solution of saddle point problems. Acta Numerica 2005; 14; 1-137. 

4. Dohrmann CR, Lehoucq RB. A primal-based penalty preconditioner for elliptic saddle point systems. SIAM Journal 
on Numerical Analysis 2006; 44(l);270-282. 

5. Tu X. A BDDC algorithm for mixed formulation of flow in porous media. Electronic Transactions on Numerical 
Analysis 2005; 20;164-179. 

6. Tu X. A BDDC algorithm for flow in porous media with a hybrid finite element discretization. Electronic 
Transactions on Numerical Analysis 2007; 26:146-160. 

7. Vassilevski PS. Multilevel block factorization preconditioners: matrix-based analysis and algorithms for solving 
finite element equations. Springer, 2008. 

8. Elman HC, Silvester DJ, Wathen AJ. Einite elements arulfa.st iterative solvers: with applications in incompressible 
fluid dynamics. Numerical Mathematics and Scientific Computation, Oxford University Press, 2005. 

9. Toselli A, Widlund OB. Domain Decomposition Methods—Algorithms and Theory, Springer Series in 
Computational Mathematics, vol. 34. Springer-Verlag, 2005. 

10. Cros JM. A preconditioner for the Schur complement domain decomposition method. Domain Decomposition 
Methods in Science and Engineering, Herrera I, Keyes DE, Widlund OB (eds.). National Autonomous University 
of Mexico (UNAM), Mexico, 2003; 373-380. 14th International Conference on Domain Decomposition Methods, 
Cocoyoc, Mexico, January 6-12, 2002. 

11. Dohrmann CR. A preconditioner for substructuring based on constrained energy minimization. SIAM J. Sci. 
Comput. 2003; 25(l):246-258. 

12. Fragakis Y, Papadrakakis M. The mosaic of high performance domain decomposition methods for structural 
mechanics: Formulation, interrelation and numerical efficiency of primal and dual methods. Computer Methods 
in Applied Mechanics and Engineering 2003; 192:3799-3830. 

13. Mandel J, Sousedlk B. BDDC and FETI-DP under minimalist assumptions. Computing 2007; 81:269-280. 

14. Sousedrk B, Mandel J. On the equivalence of primal and dual substructuring preconditioners. Electronic 
Transactions on Numerical Analysis 2008; 31:384^02. 

15. Li J, Widlund OB. BDDC algorithms for incompressible Stokes equations. SIAM Journal on Numerical Analysis 
2006; 44(6):2432-2455. 

16. Tu X, Li J. A balancing domain decomposition method by constraints for advection-diffusion problems. 
Communications in Applied Mathematics and Computational Science 2008; 3(l):25-60. 

17. Tu X. Three-level BDDC in three dimensions. SIAM Journal on Scientific Computing 2007; 29(4): 1759-1780. 

18. Mandel J, Sousedik B, Dohrmann CR. Multispace and multilevel BDDC. Computing 2008; 83(2-3):55-85. 

19. Klawonn A, Rheinbach O, Widlund OB. An analysis of a FETI-DP algorithm on irregular subdomains in the plane. 
SIAM Journal on Numerical Analysis 2008; 46(5):2484-2504. 

20. Mandel J, Sousedik B. Adaptive selection of face coarse degrees of freedom in the BDDC and the FETI-DP iterative 
substructuring methods. Computer Methods in Applied Mechanics and Engineering 2007; 196(8): 1389-1399. 

21. SIstek J, Certikova M, Burda P, Novotny J. Face-based selection of corners in 3D substructuring. Mathematics and 
Computers in Simulation 2012; 82(10):1799-181L 

22. Sousedik B, SIstek J, Mandel J. Adaptive-Multilevel BDDC and its parallel implementation. Computing 2013; 
95(12):1087-1119. 

23. Oh DS, Widlund OB, Dohrmann CR. A BDDC algorithm for Raviart-Thomas vector fields. Technical Report TR- 
951, Courant Institue of Mathematical Sciences Feb 2013. Department of Computer Science. 

24. Sousedik B. Nested BDDC for a saddle-point problem. Numerische Mathematik 2013; 125(4):761-783. 

25. Tu X. A three-level BDDC algorithm for a saddle point problem. Numerische Mathematik 2011; 119(1): 189-217. 

26. Maryska J, Rozloznik M, Tuma M. Mixed-hybrid finite element approximation of the potential fluid flow problem. 
Journal of Computational and Applied Mathematics 1995; 63:383-392. 

27. Oden J, Lee J. Dual-mixed hybrid finite element method for second-order elliptic problems. Mathematical Aspects 
of Finite Element Methods, Lecture Notes in Mathematics, vol. 606, Galligani I, Magenes E (eds.). Springer Berlin 
Heidelberg, 1977; 275-291. 

28. Maryska J, Rozloznik M, Tfima M. Schur complement systems in the mixed-hybrid finite element approximation 
of the potential fluid flow problem. SIAM Journal on Scientific Computing 2000; 22(2):704-723. 

29. Maryska J, Severyn O, Vohrallk M. Numerical simulation of fracture flow with a mixed-hybrid FEM stochastic 
discrete fracture network model. Computational Geosciences 2005; 8:217-234. 

30. Bfezina J, Hokr M. Mixed-hybrid formulation of multidimensional fracture flow. Proceedings of the 7th 
international conference on Numerical methods and applications, NMATO, Springer-Verlag: Berlin, Heidelberg, 
2011; 125-132. 

31. Bear J. Dynamics of Fluids in Porous Media. Dover, 1988. 

32. Chen Z, Huan G, Ma Y. Computational Methods for Multiphase Flows in Porous Media. SIAM, 2006. 

33. Brezzi F, Fortin M. Mixed and Hybrid Finite Element Methods. Springer-Verlag, 1991. 

34. Cowsar LC, Mandel J, Wheeler MR Balancing domain decomposition for mixed finite elements. Mathematics of 
Computation 1995; 64(211):989-1015. 


Preprint version. 



BDDC FOR POROUS MEDIA WITH COMBINED MESH DIMENSIONS 


27 


35. Pultarova I. Preconditioning of the coarse problem in the method of balanced domain decomposition by constraints. 
Math. Comput. Simulation 2012; 82(10):1788-1798. 

36. Mandel J, Dohrmann CR, Tezaur R. An algebraic theory for primal and dual substructuring methods by constraints. 
Appl. Numer. Math. 2005; 54(2): 167-193. 

37. Mandel J, Dohrmann CR. Convergence of a balancing domain decomposition by constraints and energy 
minimization. Numerical Linear Algebra with Applications 2003; 10(7):639-659. 

38. Certlkova M, Slstek J, Burda P. On selection of interface weights in domain decomposition methods. Proceedings 
of Programs and Algorithms of Numerical Mathematics 16, Dolni Maxov, Czech Republic, June 3-8, 2012. Institute 
of Mathematics AS CR, 2013; 35-44. 

39. Karypis G, Kumar V. A fast and high quality multilevel scheme for partitioning irregular graphs. SIAM Journal on 
Scientific Computing 1998; 20(l):359-392. 

40. Amestoy PR, Duff IS, L’Excellent JY. Multifrontal parallel distributed symmetric and unsymmetric solvers. 
Computer Methods in Applied Mechanics and Engineering 2000; 184:501-520. 

41. Mandel J, Sousedlk B, Slstek J. Adaptive BDDC in three dimensions. Mathematics and Computers in Simulation 
2012; 82(10):1812-1831. 


Preprint version. 



